Questioning the existence of a unique ground state structure for Si clusters 
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Density functional and quantum Monte Carlo calculations challenge the existence of a unique 
ground state structure for certain Si clusters. For Si clusters with more than a dozen atoms the 
lowest ten isomers are close in energy and for some clusters entropic effects can change the energetic 
ordering of the configurations. Isotope pure configurations with rotational symmetry and symmet- 
ric configurations containing one additional isotope are disfavored by these effects. Comparisons 
with experiment are thus difficult since a mixture of configurations is to be expected at thermal 
equilibrium. 



The determination of the structure of clusters is a dif- 
ficult task. The standard experimental techniques such 
as X-ray diffraction and NMR methods that allow to de- 
termine the atomic positions in crystals and molecules 
are not applicable to clusters . The main source of ex- 
perimental information, ion mobility measurements Q, 
gives only crude information about the overall shape of a 
cluster. The exact atomic positions of all the atoms form- 
ing the cluster remain unknown. For this reason compu- 
tational simulations provide a viable alternative to the 
experimental approach, which has been widely used for 
silicon clusters. From the theoretical point of view the 
ground state structure of a solid state system is deter- 
mined by the global minimum of the Born-Oppenheimer 
potential energy surface. Finding the global minimum 
requires global optimization algorithms. Two problems 
arise in this context. First, most global optimization al- 
gorithms give no guarantee for finding the global mini- 
mum within a finite amount of computer time. Second, 
the Born-Oppenheimer energy surface has to be calcu- 
lated with very high precision. 

Concerning the first point there is now a large amount 
of agreement between different methods for medium size 
clusters containing up to 19 atoms 0, Genetic al- 
gorithms [1, 0, Q, the big-bang method the basin 
hopping method 0, 0, OS an d the minima hopping 
method give typically similar or even identical re- 
sults. The discrepancies are rather due to different 
exchange-correlation functionals in different investiga- 
tions 

The existence of a well defined ground state structure 
is generally taken to be granted for silicon clusters. Sili- 
con clusters are, however, very different from bulk silicon 
where the second lowest configuration (a fourfold coordi- 
nated defect 14]) is 2.4 eV higher than the crystalline 
ground state. Clusters are frustrated systems, where 
most of the atoms cannot adopt their favorite fourfold co- 
ordination This can lead to small energy differences 
between different configurations. The significant devia- 
tions of the clusters bond lengths from the crystalline 



bond lengths shown in Fig. Q illustrate this frustration. 
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FIG. 1: The bond- length distribution averaged over various 
low-lying Sil7-configurations. The two vertical lines indicate 
the 1st and the 2nd nearest neighbor distances in the crystal. 



In this work we did not only search for the ground 
state configurations of silicon clusters with up to 19 
atoms, but for a large number of low energy configu- 
rations. This is possible with the dual minima hopping 
method (DMHM) |l5(, which has the property that it 
explores higher and higher energy configurations after 
having found the global minimum. Fig. [21 shows the first 
major result of our investigation, the energies of the 10 
lowest configurations of silicon clusters containing 7 to 19 
atoms. The energy difference between the global mini- 
mum and the second lowest minimum is 0.8 mHa for Sin, 
0.9 mHa for Sii3, 2.1 mHa for Sii4, 3.1 mHa for Sii7 and 
3.2 mHa for Siig. For Sii3 and Sii7 the 10 lowest con- 
figurations are in an interval of roughly 10 mHa. Since 
room temperature corresponds to ~1 mHa, entropic ef- 
fects play an important role for these clusters. 

The results of Fig. were obtained with the PBE 
functional. Even though this functional is considered to 
be among the most accurate ones, its accuracy is clearly 
insufficient to determine unambiguously the energetic or- 
dering of the configurations. For this reason we have 
performed the most accurate electronic structure calcu- 
lations that are feasible for these systems, namely quan- 
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FIG. 2: The dependence of the PBE energy interval for the 
lowest 10 configurations on the cluster size. 



turn Monte Carlo (QMC) simulations. The QMC calcu- 
lations are performed using the CHAMP code developed 
by Umrigar and Filippi. The Is, 2s and 2p electrons of 
Si are eliminated using a relativistic Hartree-Fock pseu- 
dopotential |22j |. A Slater- Jastrow type wave function is 
used as the trial wave function. The orbitals of the Slater 
determinant are taken from a DFT calculation with the 
GAMESS code using the B3LYP functional. The 
parameters of the Jastrow function describing electron- 
electron, electron-nuclear and electron-electron-nuclear 
correlations are optimized in variational Monte Carlo 
using energy minimization [24|. Diffusion Monte Carlo 
method calculates the final energies, which are presented 
in Table ffl The corresponding configurations are shown 
in Fig- 01 The QMC energies have error bars of the order 
of 1 mHa which is just enough to discriminate between 
the different energies. Even though the Monte Carlo re- 
sults change the energetic ordering of the PBE results, 
the central feature remains. Different configurations have 
energies that are nearly identical. Table [I] also shows that 
the various high quality basis sets used by different elec- 
tronic structure programs give slightly different answers 
that might change the energetic ordering. 

The new low- lying structures Sii6 a , Sii66, Sii7 a , Sim, 
Sii8a and Siig a which were found with DMHM and the 
reference structures Sii6 [ill. Sii? Siis 0] and Sii9 Q 
were already presented in [15j . The structure Sii3 was 
found by Ho et al. 6], the rotationally symmetric Sii3^ 
structure was recently proposed by Hartke [25j |. S113/ by 
Rothlisberger et al. |2£| and Sii3^ by Jeong et al. 27]. Us- 
ing DMHM ^| we found new low- lying structures Sii3 a , 
Sii36, Sii3 C and Sii3 e . From the QMC results in Tabled 
we conclude that the Sii3^ and Sii3 a configurations are 
the lowest energy structures. The new Sii3 a structure 
found with DMHM contains the stable Si6 subunit |28q . 

After having discussed the limitations of computa- 



TABLE I: The energy differences in mHa between the low 
energy geometries Sii 3 a, Sii 3b , Sii 3c , Sii 3d , Sii 3e , Sii 3 /, Sii 6 a, 
Sii6b, Sii7a, Sii7b, Sii8a, Sii9a and the reference structures 
Sii3 , Sii6, Sii7, Siis and Siig, proposed earlier as global 
minima in DFT. The Gaussian [13 calculations used the 6- 
311G(2d) basis and the DMol3 2005 Ed calculations the 
extended basis set. The CPMD |H calculations were per- 
formed with an accurate pseudopotential |21| with a 35 Ry 
plane wave cutoff and a 30 A simulation cell. 
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tional approaches in determining the total energy of sil- 
icon clusters with the necessary accuracy let us discuss 
the physical effects that can change the energetic order- 
ing. For the Sii3, Sii3 a and Sii3^ clusters we have zero 
point energies of 24.0, 25.0 and 24.5 mHa. For the Siig 
and Sii9 a we have 38.8 and 38.0 mHa. So the differences 
of the zero point energies are all of the order of mHa and 
thus not negligible, but do not change the energetic or- 
dering for the clusters we studied. In order to study the 
entropic effects we calculated the the rotational and vi- 
brational free energy based on the harmonic frequencies 
obtained from density functional (PBE) calculations [22| . 
The translational free energy does not depend on the 
configuration and was therefore not considered. If one 
compares the sum of the rotational and vibrational free 
energy for non-symmetric configurations, one typically 
finds differences of about 0.5 mHa at room temperature 
and about 1 mHa close to the melting point of the clus- 
ters |2l3. This might change the energetic ordering, but 
we did not find a case where it actually does. 

The situation is different if one compares a symmet- 
ric with a non-symmetric configuration. Silicon occurs 
in nature mainly as 28 Si or 29 Si isotope. The predomi- 
nant isotope for silicon 28 Si (abundance ~92% 31]) has 
mass 28 and no nuclear spin, the 29 Si isotope (abundance 
~5% 31]) has mass 29 and nuclear spin 1/2. When 
studying configurations with rotational symmetry, we 
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FIG. 3: Symmetrized geometries of low- lying Sii3, Sii6, Sii7, 
Sii8 and Siig isomers. 



will consider pure clusters consisting only of 28 Si atoms 
since the presence of a 29 Si atom would destroy the ro- 
tational symmetry. One can easily estimate from the 
abundancies of the isotopes that ^34% of Sii3 clusters 
will be pure clusters. For such a cluster with rotational 
symmetry, the order of the rotational subgroup enters 
into the formula for the rotational free energy. This 
leads to a weaker decrease of the free energy for sym- 
metric configurations compared to non-symmetric con- 
figurations and thus favors non-symmetric structures. In 
Fig. 2] we present the free energy curves for the structures 
Sii3 a and Sii3^ as a function of temperature with the Sii3 
free energy chosen as reference energy. The width of the 
bands for Sii3 a and Sii3d represents the statistical errors 
in the QMC energies with respect to that of structure 
Sii3. For the symmetric Sii3 a configuration the order of 
the rotational subgroup is 3, for Sii3d it is 2 and for Sii3 
it is 1. This leads to a reversal of the energetic ordering 
of the structures Sii3 and Sii3 a in the interval between 
250 and 650 K. Because of the entropic effect the Sii3 
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FIG. 4: The sum of the electronic (QMC, with errors), rota- 
tional and vibrational (including zero point) free energy con- 
tributions for Sii3a (upper band) and Sii3d (lower band) con- 
figurations as a function of temperature with the Sii3 free 
energy chosen as reference energy (solid line, set to 0). 



configuration, which is the highest at zero temperature, 
becomes the lowest at temperatures above 1000 K. At 
room temperature the Sii3 a and Sii3^ bands are sepa- 
rated by an energy gap in the range between ~1.2mHa 
and ~5.5 mHa. This corresponds to a Boltzmann weight 
in the range between 0.7% and 30%. 

These considerations are only valid for clusters consist- 
ing purely of 28 Si atoms. The presence of a 29 Si isotope 
destroys the rotational symmetry. One can estimate from 
the abundancies of the isotopes that ~24% of Sii3 clus- 
ters will contain one 29 Si isotope. If one 28 Si atom with 
nuclear spin is replaced by a 29 Si isotope which has spin 
1/2, the nuclear wave- function is a doublet and additional 
degeneracy comes from the fact that the isotope can re- 
place any of the atoms. For a non-symmetric cluster with 
N atoms the degeneracy is thus 27V. For a symmetric 
cluster that has several equivalent atoms the degeneracy 
is however reduced. In the case of the Sii3 a structure 
there are for instance only 5 non-equivalent sites, Sii3^ 
has 6 and Sii3 has 9. The nuclear entropy thus favors 
Sii3 over Sii3d by — fcTln(|) which is ~ 0.4 mHa at room 
temperature. In addition, the vibrational and rotational 
entropy contributions are slightly changed by the pres- 
ence of an isotope leading to an effect of the same order 
of magnitude. 

Up to now we have concentrated on the 10 lowest struc- 
tures. Considering higher lying configurations, the ener- 
getic spacing between configurations decreases even fur- 
ther. This can be inferred from the fact the the configu- 
rational density of states, defined as the number of con- 
figurations per energy interval, increases strongly. This 
is shown in Fig. for the Sii7 cluster. 

Several simulations have shown that the lowest energy 
structures for Si clusters with less than 20 atoms are 
non-spherical whereas larger clusters prefer to be spher- 
ical U |3 . In these studies the non-spherical to spheri- 
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FIG. 5: The configurational density of states for Sirr. The 
inset shows the ratio of the largest and smallest eigenvalues 
of the moment of inertia tensor for various low-lying Sii7 con- 
figurations. 



cal transition was obtained by considering the putative 
ground state configurations. The inset in Fig. shows 
the ratio between the largest and smallest eigenvalues of 
the moment of inertia tensor for various low- lying Sii7 
clusters. A ratio of 1 corresponds to a spherical geome- 
try, while larger values correspond to non-spherical struc- 
tures. For a given cluster size we observe in Fig. that 
the occurrence of non-spherical and spherical structures 
is independent of energy. 

In summary, we have shown that there exists a large 
number of configurations for certain silicon clusters that 
are energetically extremely close to the ground state. 
This feature was observed for Sii3 and Siig and it will pre- 
sumably be even more important for larger cluster sizes 
that were not studied in this work. As a consequence, 
entropy effects that are usually neglected can change the 
energetic ordering of the lowest configurations. Entropy 
disfavors symmetric clusters Si n in the range 13 < n < 19 
which contain in most cases no 29 Si isotope or one 29 Si 
isotope. Larger clusters will on average contain more 
than one 29 Si isotope and the symmetry related effects 
discussed above do not exist. However, for larger clusters 
the 10 lowest configurations can be expected to lie within 
an even smaller interval. The entropic effects not related 
to symmetry considerations might thus easily change the 
energy order of clusters with more than 19 atoms. Even 
if there is no reordering, different structures can be so 
close in free energy that a mixture of two or more con- 
figurations will be found at thermal equilibrium. As a 
consequence measured properties of clusters can be some 
average of the properties of several low- lying isomers. 
We thank Markus Meuwly for discussions, X. C. Zeng 
and A. Tekin for the Sii3 cluster data, the Swiss National 
Science Foundation for financial support and the staff 



of the computing center at the University of Basel for 
technical support. The work at Cornell University was 
supported by NSF grant EAR-0530301. Computational 
resources were provided by OSC, NERSC and NCSA. 



[1] B. Hartke, Angew. Chem. Int. Ed., 41, 1468 (2002). 

[2] R. R. Hudgins, M. Imai, M. F. Jarrold and P. Dugourd, 

J. Chem. Phys. Ill, 7865 (1999). 
[3] X. Zhu and X. C. Zeng, J. Chem. Phys. 118, 3558 (2003). 
[4] X. L. Zhu et al, J. Chem. Phys. 120, 8985 (2004). 
[5] B. Hartke, J. Phys. Chem. 97, 9973 (1993). 
[6] K.-M. Ho et al, Nature 392, 582 (1998). 
[7] I.Rata et al, Phys. Rev. Lett. 85, 546549 (2000). 
[8] K. A. Jackson et al, Phys. Rev. Lett. 93, 013401 (2004). 
[9] Z. Li, H. Scheraga, Proc. Natl. Acad. Sci. USA 84, 6611 

(1987). 

[10] J. Doye and D. Wales, Phys. Rev. Lett. 80, 1357 (1998). 
[11] S. Yoo, X. C. Zeng, Angew. Chem. Int. Ed. 44, 1491 
(2005). 

[12] S. Goedecker, J. Chem. Phys. 120 9911 (2004). 
[13] S. Yoo and X. C. Zeng, J. Chem. Phys. 123, 164303 
(2005). 

[14] S. Goedecker, T. Deutsch and L. Billard, Phys. Rev. Lett. 
88, 235501 (2002) 

[15] S. Goedecker, W. Hellmann and T. Lenosky, Phys. Rev. 
Lett. 95, 055501 (2005). 

[16] J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett. 
77, 3865 (1996). 

[17] Gaussian 03, Revision B.01, M. J. Frisch, G. W. Trucks, 
H. B. Schlegel, G. E. Scuseria et al Gaussian, Inc., Pitts- 
burgh PA, 2003. 

[18] B. Delley, J. Chem. Phys. 92, 508 (1990). 

[19] B. Delley, J. Chem. Phys. 113, 7756 (2000). 

[20] CPMD Version 3.3: developed by J. Hutter, A. Alavi, 
T. Deutsch, M. Bernasconi, S. Goedecker, D. Marx, M. 
Tuckerman and M. Parrinello, Max-Planck-Institut fur 
Festkorperforschung and IBM Zurich Research Labora- 
tory (1995-1999). 

[21] S. Goedecker, M. Teter, J. Hutter, Phys. Rev. B 54, 1703 
(1996). 

[22] J. R. Trail and R. J. Needs, J. Chem. Phys. 122, 174109 

(2005). 

[23] http://www.msg. ameslab.gov/GAMESS/ 
[24] C. J. Umrigar and C. Filippi. Phys. Rev. Lett. 94, 150201 
(2005). 

[25] A. Tekin, B. Hartke, Phys. Chem. Chem. Phys. 6, 503 

(2004) . 

[26] U. Rothlisberger et al, J. Chem. Phys. 96, 1248 (1992). 
[27] J. Jeong et al, J. Phys.: Condens. Matter 10, 5851 
(1998). 

[28] M. Jarrold and E. Bower, J. Phys. Chem 92, 5702 (1988). 

[29] The formulas for the rotational and vibrational free en- 
ergy are given in standard textbooks, see for example 
F.Jensen, Introduction to Computational chemistry (Wi- 
ley, New York, 1999). 

[30] F. Baletto and R. Ferrando, Rev. Mod. Phys. 77, 371 

(2005) . 

[31] T. B. Coplen et al, Pure Appl.Chem. 74, 1987 (2002). 



